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ABSTRACT 

Massive stars can be efficiently ejected from their birth clusters through encounters with other massive stars. 

We study how the dynamical ejection fraction of O star systems varies with the masses of very young star 
clusters, Med, by means of direct A^-body calculations. We include diverse initial conditions by varying the 
half-mass radius, initial mass-segregation, initial binary fraction and orbital parameters of the massive binaries. 

The results show robustly that the ejection fraction of O star systems exhibits a maximum at a cluster mass of 
10^-® Mq for all models, even though the number of the ejected systems increases with cluster mass. We show 
that lower mass clusters (Med ~ 400 Mq) are the dominant sources for populating the Galactic field with 
O stars by dynamical ejections, considering the mass function of embedded clusters. About 15 per cent (up 
to ~38 per cent, depending on the cluster models) of O stars of which a significant fraction are binaries, and 
which would have formed in a rjIO Myr epoch of star formation in a distribution of embedded clusters, will 
be dynamically ejected to the field. Individual clusters may eject 100 per cent of their original O star content. 

A large fraction of such O stars have velocities up to only 10 km s“^. Synthesising a young star cluster mass 
function it follows, given the stellar-dynamical results presented here, that the observed fractions of field and 
runaway O stars, and the binary fractions among them can be well understood theoretically if all O stars form 
in embedded clusters. 

Keywords: methods: numerical - stars: kinematics and dynamics - stars: massive- open clusters and associa¬ 
tions: general - galaxies: star clusters: general 


1. INTRODUCTION 

The majority of O stars are found in young star clusters 
or OB associations. It has been suggested that a small frac¬ 
tion of massive stars may form in isolation as a result of 
stochastic sampling of the stellar and cluster mass function 
(see de Wit et al. 2005; Parker & Goodwin 2007; Oey & Lamb 
2012). However, in the case of the Galactic field O stars 
within 2-3 kpc from the Sun, essentially all of them can be 
traced back to young star clusters or associations (Schilbach 
& Roser 2008), suggesting the field O stars also formed in 
clusters (Gvaramadze et al. 2012). Some observational stud¬ 
ies on O stars in the Large and Small Magellanic Clouds have 
claimed there to be evidence for isolated massive star forma¬ 
tion (Bressert et al. 2012; Oey et al. 2013), such as e.g., ap¬ 
parent isolation, absence of a bow-shock, and circular shape 
of an H II region. These claims can not, however, be conclu¬ 
sive since they can be explained with slowly moving runaways 
or former members of clusters which have already dissolved 
(Weidner et al. 2007). Eor example, only 30^0 per cent of 
the Galactic OB runaways are found to have a bow-shock (van 
Buren et al. 1995; Huthoff & Kaper 2002), and, furthermore, 
bow-shock formation depends on the physical condition of the 
ambient medium through which a runaway travels (e.g., tem¬ 
perature and density of the ambient medium, Huthoff & Kaper 
2002). Mackey et al. (2013) showed that runaways can have 
a circular H II region in the projected Ha emission. Eurther- 
more, in the case of the 30 Dor region Bressert et al. consid¬ 
ered, there are many young star clusters in the region. If the 
30 Dor region is modelled by a population of clusters, how 
many O stars would be ejected or removed from their clusters 
in 10 Myr? This has not been investigated. Any conclusion on 
the issue of isolated O-star formation are, for the time being, 
thus premature. Einally, measuring stellar proper motions in 
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the Large or Small Magellanic Clouds is unfeasible such that 
claims for O stars formed in isolation become next to impos¬ 
sible to confirm or reject. With Gaia, however, it may become 
possible to put constraints on this if proper motions of mas¬ 
sive stars can be measured to a precision of a few km in 
the Magellanic Clouds. Given the above results on the well 
studied ensemble of O stars within 2-3 kpc around the Sun 
it appears physically more plausible that all O stars form in 
embedded clusters. 

There are two mechanisms to expel O stars from a star clus¬ 
ter. The one is the binary supernova scenario (Blaauw 1961; 
Portegies Zwart 2000; Eldridge et al. 2011) in which a star 
in a binary obtains a high velocity after the supernova ex¬ 
plosion of the initially more massive companion star. The 
other is dynamical ejection through strong close encounters 
with other massive stars in the cluster core (Poveda et al. 
1967; Leonai'd & Duncan 1990; Eujii & Portegies Zwart 2011; 
Banerjee et al. 2012; Perets & Subr 2012). Of the candidate 
isolated O stars catalogued by Schilbach & Roser (2008), 8- 
9 per cent could not be traced back to any still existing star 
clusters or associations, but are consistent with originating 
in a cluster which could not be found due to several effects 
such as incompleteness of their cluster sample (Schilbach & 
Roser 2008), dissolution of the cluster (Gvaramadze et al. 
2012), and having experienced the two-step-ejection process 
proposed by Pflamm-Altenburg & Kroupa (2010). The ma¬ 
jority of the stars which can be traced back to clusters appear 
to have been ejected when their parental clusters were very 
young (Schilbach & Roser 2008). This may suggest that dy¬ 
namical ejections through few-body interactions in the heavy- 
star-rich cluster core mainly populate the field O star popula¬ 
tion (Clarke & Pringle 1992). 

There have been many studies on dynamically ejected mas¬ 
sive stars from young star clusters using numerical integra¬ 
tions. In most cases the studies performed few-body scat- 
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tering experiments focused on runaways, which are ejected 
with a velocity higher than rs 30^0 km (Gualandris et al. 
2004; Pflamm-Altenburg & Kroupa 2006; Gvaramadze et al. 
2009; Gvaramadze & Gualandris 2011). Full 7V-body calcu¬ 
lations for a particular cluster mass have also been presented 
(Banerjee, Kroupa, & Oh 2012; Perets & Subr 2012). Fuji! 
& Portegies Zwart (2011) performed A^-body calculations for 
very dense star clusters composed of single stars with three 
different cluster masses (> 6000 Mq) and showed that the 
clusters with the lowest mass in their models, 6000 Mq, pro¬ 
duce runaways most efficiently. They argued that the one bul¬ 
lying binary formed during core collapse of the cluster is re¬ 
sponsible for the production of runaways and that the number 
of runaways is almost independent of cluster mass. The run¬ 
away fraction (number of runaway stars over total number of 
stars still present) therefore decreases with cluster mass as the 
more massive clusters have more stars. Note that the lower- 
mass limit for stars in their models is 1 Mq. Thus, if a cluster 
is populated with a full range of stellar masses from the hydro¬ 
gen burning limit upwards, in order to have the same amount 
of massive stars, the cluster would be twice as massive than 
they claimed. The relation between O star ejection fraction 
and cluster mass has not been studied in detail to date. 

Here we study for the first time how the dynamical ejection 
fraction of O stars changes with the masses of very young star 
clusters using direct A^-body calculations with a broad range 
of cluster masses, from low-mass clusters containing one or 
two O stars to massive clusters containing a few hundred O 
stars, and including primordial binaries. Such an investigation 
based on a very large library of young cluster models has not 
been possible to date because the knowledge on the proper¬ 
ties of initial stellar population, their birth configurations and 
last but not least the computational algorithms enabling such 
a CPU-intensive and massive numerical challenge have not 
been in place. Thus, for example, how to initialize initially 
a mass segregated cluster in dynamical equilibrium is only 
known since the work of Subr et al. (2008) and Baumgardt 
et al. (2008) and the statistical distribution functions of mas¬ 
sive star binaries are becoming observationally constrained 
only now (e.g., Sana et al. 2012,2013; Kiminki & Kobulnicky 
2012; Kobulnicky et al. 2014; Aldoretta et al. 2015). The N- 
body calculations, made possible by the relentless code and 
algorithmic advances by Sverre Aarseth and Seppo Mikkola, 
used in this study are briefly described in Section 2. The re¬ 
sults are presented in Sections 3-6, and the discussion and 
summary follow in Sections 7 and 8, respectively. 

2. A-BODY MODELS 

The data used in this study are based on the theoretical 
young star cluster library of model sequences over cluster 
mass computed by Oh & Kroupa (2012) with direct A^-body 
calculations using the NBODY6' code (Aarseth 2003) and 
various initial conditions. Among the model sequences in 
the library, we adopt four sequences here and, further, extend 
them to higher cluster masses. Additionally we perform cal¬ 
culations for two more sets of initial conditions for this study. 
Table 1 lists properties of the model sequences studied in this 
paper. The initial setup of the models is described below. 

The four model sequences that we adopt from the library 
for this study are comprised of binary-rich clusters with initial 
mass segregation (MS30P) and without initial mass segrega- 

* The code can be freely downloaded from 
http://www.ast.cam.ac.uk/'-sverre/web/pages/nbody.htm. 


Table 1 

List of A^-body model sequences studied here. 


Name 

rh (pc) 

IMS 

IBF 

IPD 

Pairing method 

MS30P 

0.3 

Yes 

1 

Kroupa95 

OP 

NMS30P 

0.3 

No 

1 

Kroupa95 

OP 

MS30P.SP 

0.3 

Yes 

1 

Sana et al.l2 

OP 

MS3UQ.SP 

0.3 

Yes 

1 

Sana et al.l2 

uniform q dist. 

MS3S 

0.3 

Yes 

0 



MS80P 

0.8 

Yes 

1 

Kroupa95 

OP 

MSI OP 

0.1 

Yes 

1 

Kroupa95 

OP 


Note. — Names of model sequences and initial half-mass radii (rh) 
are listed in columns 1 and 2, respectively. Columns 3-6 present initial 
mass segregation (IMS), initial binary fraction (IBF), and initial period 
distribution (IPD) and pairing method used for massive binaries (primary 
mass > 5 Mq). 


tion (NMS30P), single-star clusters (MS3S) with an initial 
half-mass radius, rh, of 0.3 pc, and initially mass-segregated, 
binary-rich clusters with an initial rh of 0.8 pc (MS80P). 

The cluster masses, Med, in the library range from 10 to 
10^-^ Mq (to 10^ Mq for the MS30P model sequence). For 
extending the range of cluster masses, we additionally carried 
out calculations for clusters with Med ~ 10^ and lO"*’® Mq 
(Med ~ 10^'^ Mq for the MS30P model sequence). Indi¬ 
vidual stellar masses are randomly drawn from the canonical 
initial mass function (IMF), which is a two-part power-law 
(Kroupa 2001; Kroupa et al. 2013), 


^(m) = k 


{^)y, _ 0.08 < m/Me < 0.50, 

(sl) ' '(io)’'''’ 0.50 < m/MQ < 


( 1 ) 

More details of our stellar mass sampling can be found in 
Oh & Kroupa (2012). The mass of the most massive star 
in the cluster, mmax, is chosen from the maximum-stellar- 
mass - cluster-mass relation (Weidner & Kroupa 2004; Wei- 
dner et al. 2010, 2013a, 2014). Thus only clusters with 
Med > 10^'® Mq initially have O-type stars in our models. 
Pflamm-Altenburg et al. (2007) provide a fitting formula for 
Wniax(Vfed)- This procedure was adopted rather than opti¬ 
mal sampling introduced in Kroupa et al. (2013) because the 
theoretical young cluster library of Oh & Kroupa (2012) did 
not have this sampling method to disposal. 

Binary-rich models have an initial binary fraction of 100 
per cent, i.e., all stars are in a binary system at f = 0 Myr. 
This is motivated by the angular momentum problem of star 
formation since the angular momentum of a collapsing cloud 
core can be distributed efficiently into two stars, and also by 
the strong empirical evidence that star formation universally 
prefers the binary mode (Duchene et al. 2007; Goodwin & 
Kroupa 2005; Goodwin et al. 2007; Marks & Kroupa 2012; 
Leigh et al. 2015; Sana et al. 2014). The Kroupa (1995b) 
period distribution function with minimum and maximum pe¬ 
riods of 10 and 10® "^^ days, respectively, is adopted for all 
binaries in the library. The Kroupa period distribution. 


.fp = 2.5 


logio P - 1 

45+(logioU-l)2’ 


( 2 ) 


where the period (P) is in days, has been derived iteratively 
using binary-star data from the Galactic field population and 
from star forming regions. With this distribution, about 2.75 
per cent of all O star binaries have a period shorter than 
100 days (Figure 1). Although it was derived from binaries 
with primary star mass < 1 Mq, the same distribution is used 
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for more massive binaries in the library since no well con¬ 
strained period distribution of massive binary population was 
available when the 7V-body calculations were begun. 

However, recently, Sana et al. (2012) derived an 
intrinsic period distribution of O stars of the form 
fp oc (log^o over the period range of 0.15 < 

log]^Q(P/days) < 3.5, from spectroscopic observations of 
O-star populations in nearby young open star clusters. The 
Sana et al. (2012) period distribution leads to « 49 per cent 
of O stars being binaries with a period shorter than 100 d 
(Figure 1). Here we additionally perform a set of A^-body 
integrations with massive binaries (m > 5 Mq) having the 
Sana et al. (2012) period distribution, MS30P_SP, but with the 
maximum period extended to the value log]^Q(P/days) « 6.7 
at which the cumulative binary fraction becomes unity. We 
thus adopt the following period distribution function for the 
massive binaries, 

/P = 0.23x (logloP)-°•"^ (3) 

for which fp d log^o P = 1. The other initial conditions 
are the same as those of MS30P. See also Kiminki & Kob- 
ulnicky (2012), Kobulnicky et al. (2014) and Aldoretta et al. 
(2015) for slightly different period distributions of O-star bi¬ 
naries derived from observations. 

Stars with a mass m < 5 Mq are randomly paired after 
choosing the two masses randomly from the IMF. Again, this 
mass-ratio distribution is derived from the Galactic field pop¬ 
ulation and from star forming regions (Kroupa 1995a,b). Stars 
with m > 5 Mq are paired to imprint observations that the 
massive binaries favour massive companions (Pinsonneault & 
Stanek 2006; Kobulnicky & Fryer 2007; Kiminki & Kobul¬ 
nicky 2012; Sana et al. 2012). To achieve this we first order 
all stars more massive than 5 Mq in to an array of decreas¬ 
ing mass. The most massive star is then paired with the next 
massive star, the third most massive star is paired with the 
fourth most massive star, and so on. We refer to this pro¬ 
cedure as ordered pairing (OP). Whilst with this method it 
is easy to have a massive primary being paired with a mas¬ 
sive companion while preserving the IMF, it gives a signifi¬ 
cant bias toward high mass ratios (<? « 1 where q = 1712 !mi 
and TOi > m 2 ) which is not supported by recent observations 
(Kiminki & Kobulnicky 2012; Sana et al. 2012; Kobulnicky 
et al. 2014). The recent observations show the mass-ratio dis¬ 
tribution of O-type binaries to be rather uniform, 

fiq) « g" 

with —0.1 < 77 < 0.1 and 0.1 < q < 1.0. We add a model se¬ 
quence (MS3UQ_SP) with a uniform mass-ratio distribution, 
?7 = 0, for massive binaries (mi > 5 Mq), being otherwise 
the same as the MS30P_SP model sequence. First, a mass- 
ratio value for a primary is derived from a uniform probability 
function and then the secondary is chosen from the rest of the 
members which gives the closest q value for the derived one. 
This is important to preserve the IMF; It is important to first 
generate a full list of stars whose combined mass corresponds 
to the mass of the cluster and only then to create the binaries 
using this list only. The thermal distribution is used for the 
initial eccentricity distribution, i.e., /e(e) = 2e, where e is 
the eccentricity of a binary (see Kroupa 2008). 

Each model sequence has one initial half-mass radius re¬ 
gardless of the cluster mass since by observation there is no 
significant relation between the size and mass of a cluster 
(Bastian et al. 2005; Scheepmaker et al. 2007) or at most a 



Figure 1. Cumulative binary fraction for logjQ P < logj^g P' among all O 
stars. Solid and dashed lines are calculated from the Kroupa (1995b, Equa¬ 
tion (2)) and Sana et al. (2012, Equation (3)) period distributions, respec¬ 
tively. Dotted line is a uniform period distribution with a period range of 
0.5 < log]^g(P/days) < 4 used in Banerjee et al. (2012) for O star bina¬ 
ries. For example, 49 per cent of all O stars are binaries with period shorter 
than 100 d according to the Sana et al. distribution. 

only very weak one (Larsen 2004; Marks & Kroupa 2012). 
The cluster-mass-initial half-mass radius relation from Marks 
& Kroupa (2012), 

(fh/pc) = 0.1 X (Mecl/M 0 )°-^^, 

gives initial half-mass radii of « 0.21 to 0.38 pc for Med = 
102-5 _io 4-5 Q^jj. initial half-mass radius (0.3 pc) 
lies close to the median value and is well consistent with the 
typical initial half-mass radius (0.1-0.3 pc) inferred for the 
Milky Way clusters (Marks & Kroupa 2011). Thus, our re¬ 
sults may be able to represent reality. Here we include two 
more model sequences with different half-mass radii of 0.1 
and 0.8 pc, MSlOP and MS80P, respectively. This allows 
us to estimate how the ejection fraction changes with cluster 
mass when the cluster radius varies with cluster mass. While 
the MS80P model sequence is adopted from the library of 
Oh & Kroupa (2012), the MSlOP model sequence is com¬ 
puted for this work. From the choices of half-mass radii, the 
initial central densities of our models range from « 300 to 
3.3 X lO'^ Mqpc~^ for rh = 0.8 pc clusters, from 6.2 x 10^ 
to 6.2 X 10^ Mqpc~^ for rh = 0.3 pc clusters, and from 
Rs 1.7xlO®M0pc“^toft! 1.7xlO^M0pc“^forrh = 0.1 pc 
clusters. 

For all models, initial positions and velocities of single stars 
or of centre-of-masses of binaries in the clusters are gener¬ 
ated following the Plummer density profile under the assump¬ 
tion that the clusters are in virial equilibrium. Initially mass- 
segregated clusters are produced being fully mass-segregated 
with the method developed in Baumgardt, De Marchi, & 
Kroupa (2008) in which more massive systems are more 
bound to the cluster. 

The A^-body code regularises the motions of stars in com¬ 
pact configurations and incorporates a near-neighbour force 
evaluation scheme to guarantee the high accuracy of the solu¬ 
tions of the equations of motions. Stellar evolution is taken 
into account by using the stellar evolution library (Hurley 
et al. 2000, 2002) implemented in the code. All clusters are 
evolved up to 3 Myr, i.e., to a time before the first supernova 
occurs, in order to study dynamical ejections only. Moreover 
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the dynamical ejections of O stars are most efficient at an early 
age of the cluster (Oh et al. in preparation and see Section 5), 
thus computations to an older cluster age would not change 
our conclusion on the dynamical ejections. 

A standard solar-neighbourhood tidal field is adopted. The 
clusters are highly tidally underfilling such that the tidal field 
of the hosting galaxy plays no significant role for the re¬ 
sults presented here. Our cluster models are initially gas- 
free, i.e., without a background gas-potential, and thus the 
gas-expulsion phase is not included here. While 100 realisa¬ 
tions with different random seed numbers are performed for 
the clusters with Med < 10^'® Mq, 10 and 4 realisations 
are carried out for Med = 10^ and 10^'^ Mq clusters, re¬ 
spectively, as computational costs increase with the number 
of stars, i.e., cluster mass, and stochastic effects are reduced 
with a larger number of stars. 

3. EJECTION FRACTION OF O STARS AS A FUNCTION OF ATed 

Throughout this paper, an O star refers to a star with a mass 
of > 17.5 Mq at 3 Myr and an O-star system is either a single 
O star or a binary system having at least one O star compo¬ 
nent. We regard a system as being ejected when its distance 
from the cluster centre is greater than 5 x of the cluster 
at 3 Myr, and its velocity is larger than the escape velocity at 
that radius. Such a short distance criterion for ejections is an 
appropriate physical condition to quantify the real ejected O 
star fraction, in order to account for the O stars which have 
been dynamically ejected but have not travelled far from the 
cluster centre by 3 Myr. O stars can be ejected with a velocity 
of a few km from the lower-mass clusters, e.g., the cen¬ 
tral escape velocity of a cluster with Med = 1000 Mq and 
rh = 0.3 pc is « 6 km s“^. 

From the computational data the ejection fraction of O star 
systems, ’^°™fei,o, is defined as the number ratio of the ejected 
O star systems, to all O star systems that are still 

present in a calculation, Nq, i.e., 

com AT 

com £ . _ 

lej,0 com^Q 

Only clusters with Mgd > 10^'® Mq are shown in the tables 
because clusters with Med < 100 Mq do not have an O star 
in our models. Since values of ‘^°™/ej,o from individual clus¬ 
ters are diverse, especially for low-mass clusters, we mainly 
deal with the averaged value of ™™/ej,o. (““/ej,o). for the 
same cluster mass in a certain initial conditions set. 

The average ejection fractions of O star systems for all 
models are shown in Figure 2. Table Al lists averaged proper¬ 
ties of the models, such as half-mass radius, number of O star 
systems, number of ejected O star systems, O star ejection 
fraction, and O star binary fraction in the clusters and among 
the ejected O star systems, at 3 Myr. All models show a trend 
that ('^°™/ej,o) increases with a cluster mass up to 10^ ® Mq 
and then decreases for a higher cluster mass. The average val¬ 
ues of the ejection fraction with the same cluster mass vary 
from model to model. The peak fractions are « 25 per cent for 
the binary-rich cluster models with rh = 0.3 pc with a weak 
dependency on initial mass segregation or on the initial pe¬ 
riod distribution of massive binaries, although the ('^°™/ej,o) 
values of the other cluster masses change with the initial con¬ 
ditions. Note that individual clusters in the mass range 10^’® 
to 10^’® Mq may eject all their O star systems! 

The ejection fractions from the initially not mass- 
segregated models (NMS30P) exhibit similar results to the 


MS30P models. The difference between the two models is 
only seen at the lowest cluster mass, 10^'® Mq. This can be 
understood because dynamical mass-segregation is inefficient 
for the lowest mass clusters due to their low density. Initial 
mass segregation therefore does not have a significant impact 
on massive-star ejections. 

In the models using the Sana et al. (2012) period distribu¬ 
tion (MS30P_SP and MS3UQ_SP), the same trend of the O 
star ejection fraction is present as above, but the (‘^°™/ej,o) 
values show a much broader peak than those in the model 
using the Kroupa (1995b) period distribution (MS30P). The 
(comjVej^o) values of the three models are similar (see Ta¬ 
ble Al). However, for massive clusters (Med > 10^'^ Mq), 
MS30P_SP clusters not only have a smaller number of O star 
systems near 3 Myr due to a higher fraction of binaries and 
mergers, but they also eject more O star systems. The latter 
suggests that short period massive binaries are dynamically 
active and boost the ejection of O stars, especially for mas¬ 
sive clusters. 

The numbers of O-star systems are slightly higher in the 
MS3UQ_SP sequence than in the MS30P_SP sequence (see 
Table Al) because O-star binaries in the MS30P_SP mod¬ 
els mostly have an O-star companion while companions of 
O -star binaries in the MS3UQ_SP models are chosen from a 
broader range of stellar masses. For example, when the binary 
fraction is unity and both sequences have the same number of 
individual O stars, for the MS30P_SP sequence the number 
of O-star systems would be almost half of the number of in¬ 
dividual O stars while it would larger than half in the case 
of the MS3UQ_SP sequence. However, the difference in the 
mass-ratio distributions of the two model sequences does not 
seem to affect the trend of the O star ejection fraction as a 
function of cluster mass having a peak at Med = 10^'^ Mq. 
The difference between the results of the MS30P_SP and 
MS3UQ_SP models is only marginal. Therefore, massive star 
ejections are not sensitive to the form of the mass-ratio dis¬ 
tribution as long as they are paired with other massive stars, 
i.e., (j > 0.1. 

For single-star clusters (MS3S), we find the same cluster 
mass at which the peak of the O star ejection fraction oc¬ 
curs although the ejection fractions are about half those of 
the binary-rich model sequences. The peak appears at the 
same cluster mass in the sequences with different initial clus¬ 
ter sizes. In the case of MS80P (rh = 0.8 pc), the ejection 
fraction is much smaller than in other model sequences with 
rh = 0.3 pc. But (™™/ej,o) reaches « 50 per cent in the 
MSI OP model sequence (initial rh = 0.1 pc) at the cluster 
mass of the peak. The dependency of the ejection fraction on 
cluster mass is thus independent of the initial cluster radius, 
but more compact cluster do have higher ejection fraction of 
massive stars. 

The peak near 10^'® Mq of the ejection fraction can be 
understood as a consequence of the growth of the num¬ 
ber of O stars with cluster mass surpassing the number of 
ejected O stars as the cluster mass increases and the poten¬ 
tial smoothens. Figure 3 shows the average number of ejected 
O stars at each cluster mass. 

This implies that a-few-Myr old clusters with a mass at 
which the O-star ejection fraction peaks would present O 
star populations most deviating from the number of O stars 
given by the canonical IMF. For example the estimated stel¬ 
lar mass in the inner 2 pc of the Orion Nebula Cluster (ONC) 
is Rs 1800 Mq (Hillenbrand & Hartmann 1998) and proba- 
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logio (j^ec./Mo) 


Figure 2. Ejection fraction of O-star systems as a function of cluster mass at 3 Myr. Red big circles are the average O star ejection fraction for each cluster 
mass and open circles are the values of individual clusters. The (red) solid lines are drawn by connecting the red points to guide the dependency of the ejection 
fraction on the cluster mass. Grey vertical bars indicate where the central 68 per cent of the data points lie (i.e., the points between 16 (lower) and 84 (upper) 
percentiles) for 10^ and 10^'^ Mq cluster models (for the MSlOP sequence 10^'^ Mq clusters are included). Because more than 85 per cent of clusters with 
Meci = 10^'^ Mq in each sequence except for clusters in the MSlOP sequence do not eject any O star, the grey bar is not plotted for this cluster mass with 
the exception for the MSlOP sequence. For more than 50 per cent of the 10^ Mq clusters of all sequences, except for the MSlOP sequence for which only 
30 per cent of the cases, the O star ejection fraction is 0. Due to their small number of realisations and small spread of the ejection fraction, the grey bar is not 
necessary for clusters with Med ^ 10"^ Mq. Blue triangles are the average O star ejection fraction using 10 pc for the ejection criterion without applying a 
velocity criterion. The data for the 10^ Mq cluster model taken from the calculations by Banerjee et al. (2012) are marked with open stars in the MS80P model 
sequence panel. 
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bly « 2700 Mq taking account of a probable binary fraction 
of 50 per cent (Kroupa 2000). With these cluster masses, 
about 6-10 stars are expected to have a mass larger than 
18 Mq from the canonical IMF. However, only three O-type 
(or m > 17.5 Mq) stars are found within « 2 pc of the 
Trapezium stars (Hillenbrand 1997). This discrepancy may 
be due to the stochastic sampling of the IMF. But the proba¬ 
bility for this is less than 5 per cent when stellar masses are 
randomly drawn from the IMF with a stellar upper mass limit 
of 150 Mq (but about 10 per cent with TOmax = 50 Mq) 
such that it is more likely that the ONC has ejected a signif¬ 
icant fraction of its O star content. Furthermore, the initial 
mass of the cluster may have been higher than observed at the 
present day as the cluster could have lost a significant frac¬ 
tion of its members during the residual gas expulsion phase 
(Kroupa et al. 2001). Thus the ONC formed with a mass close 
to the peak mass we found in this study and could have lost a 
significant fraction of its massive stars through the dynamical 
ejection process. This may explain why the ONC’s massive 
star population deviates from the canonical IMF (Pflamm- 
Altenburg & Kroupa 2006). Examples of such events would 
be two O-type single runaway stars, AE Aur and p Col, and 
a massive binary system l Ori which may have been ejected 
from the ONC through a binary-binary interaction (Gies & 
Bolton 1986; Hoogerwerf et al. 2001; Gualandris et al. 2004). 
Euture proper motion surveys, such as with the Gaia mission, 
may provide further constraints on this question, because the 
putative ejected O stars from the ONC would have to become 
evident. 

Very massive clusters may have small ejection fractions 
of O stars but their very upper stellar mass end could lack 
stars because the ejection fraction of stars more massive than 
100 Mq can be from « 10 to « 100 per cent (Banerjee et al. 
2012). Eurthermore, they can eject not only very massive sin¬ 
gle stars but also very massive binaries (Oh et al. 2014). This 
bias through the most massive stars being ejected from star- 
burst clusters and the observed IME of star-burst clusters be¬ 
ing close to the canonical value (e.g., the R136 cluster in the 
Large Magellanic Cloud, Massey & Hunter 1998) may im¬ 
ply the IME to be top-heavy in star-burst clusters (Banerjee & 
Kroupa 2012). There is some evidence for a possibly system¬ 
atically varying IME above a few Mq : A number of indepen¬ 
dent arguments have shown the IME to become increasingly 
top-heavy with increasing star-formation rate density on a pc 
scale (Marks et al. 2012), with a dependency on the metallic- 
ity. This metallicity dependence is such that metal-rich envi¬ 
ronments tend to produce a steeper (i.e. top-light) IME slope. 
Evidence for this may have emerged in M31 data (Weisz et al. 
2015). The models calculated here are, however, in the invari¬ 
ant regime. 

A theoretically expected relation between the ejection frac¬ 
tion and the cluster mass can be found. Eirst, the number of O 
stars in a cluster, A"o,imf, which is a function of TOniax(A/eci) 
and Meci, is 


.^O.iMF as a function of Med from Equation (4). Since we 
adopt the mmax-Med relation. Equation (4) is not a simple 
linear function of cluster mass especially for low-mass clus¬ 
ters in which mmax is highly dependent on the cluster mass. 
.iVo.iMF, however, becomes almost linearly proportional to 
Med for massive clusters where mmax is saturated (see the 
upper panel of Eigure 4). Note that the solid line in the upper 
panel of Eigure 4 and the numbers in Table Al are slightly dif¬ 
ferent because Equation (4) gives the number of all individual 
O stars from the IME, whereas (‘^°™A^o) listed in Table Al 
is the number of O star systems at 3 Myr. The difference 
between the number of O stars from the IME, A^o.imf, and 
from our A^-body calculations, (“™7Vo), at the cluster mass 
of 10^'® Mq in Eigure 4 arises due to our choice of forcing a 
cluster to have at least one mmax-star initially. This would not 
have been required if optimal sampling (Kroupa et al. 2013) 
had been available when the theoretical young star cluster li¬ 
brary was computed pre-2013. 

Our stellar masses are randomly drawn from the IME with 
the upper mass limit given by the mmax-Med relation. There¬ 
fore, our sampling would produce a slightly lower maximum 
stellar mass in a cluster but a similar number of O stars com¬ 
pared to what random sampling without the mmax-Med rela¬ 
tion would give. The grey line in the upper panel of Eigure 4 is 
the expected number of O stars from the IME with a universal 
upper stellar-mass-limit of 150 Mq for all cluster masses. The 
figure indicates that a higher number of O stars is expected in 
the IME with a universal and constant mmax = 150 Mq than 
with the mmax-Med relation at lower cluster mass (Med ^ 
a few 1000 Mq, the black solid line in the same figure). But 
the values of in our A^-body models are close to the 

grey line because of adding one mmax(Med)-star in a clus¬ 
ter of mass Med by default. So whether adopting the relation 
or not would have little effect on the number of O stars in 
our models. The proper comparison for outcomes of the two 
assumptions (mmax = 150 Mq or the mmax-Med relation) 
can be done by performing a large set of Monte Carlo calcu¬ 
lations. However, it is beyond the objectives of this study. In 
any case, the mmax-Med relation is observationally well sup¬ 
ported for young star clusters and star-forming regions in the 
Milky Way given the most recent data (Kirk & Myers 2011; 
Weidner et al. 2013a). 

Eitting the number of O-star systems present, regardless of 
whether the systems are in the star cluster or ejected, in N- 
body calculations at 3 Myr to a simple power law gives. 


r^No) oc 


' ^0.97±0.01 

^0.96±0.01 

"-'ed , 1 

^l.OliO.Ol^ 

»/foo±o.oi’ 

-''^ed > 


MS30P, 

NMS30P, 

MS30P_SP, 

MS3UQ.SP, 

MS3S, 

MS80P, 

MSIOP. 


(5) 


^TTlmax (-^^ecl) 

lVo,iMF= / ^(m)dm, 

J17.5 Mq 

3.390 X 10-4 - O.Oldm-i^f , ^ 

« - max 

0.100- 0.062mmax 

by using the IMP (Equation (1)). The normalisation con¬ 
stant, k, in Equation (1) is calculated from Med = 
tnC(''7t)dm. The top panel of Eigure 4 presents 


The numbers of O star systems in all model sequences are 
more or less linearly proportional to the cluster mass as ex¬ 
pected from the IMP despite various binary fractions. 

Next, the relation between the number of ejected O stars 
and cluster mass is needed. According to Verbunt (2003), the 
binary - single star scattering encounter rate, F, is 


^cl 


( 6 ) 
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where po is the central mass density, Uci the velocity disper¬ 
sion of the cluster, Tc the core radius, and a the semi-major 
axis of the binary. The semi-major axis can be replaced with 
the system mass, mgys, and the orbital velocity, Vorb, of the 
binary by using Equation (8.142) in Kroupa (2008) as, 

a oc msysV~l. (7) 

Considering binaries with a circular velocity similar to the ve¬ 
locity dispersion of the cluster are dynamically the most active 
(i.e., Vorb = CTci) and ad oc y/poVc, Equation ( 6 ) becomes 

r oc Po^nisys oc M°^fmsys, ( 8 ) 

since all our cluster models in a given sequence have the same 
initial half-mass radius the initial central density is only pro¬ 
portional to the cluster mass. A binary ejecting an O star af¬ 
ter a strong interaction with a single star is likely a binary 
consisting of O stars since the least massive star is generally 
ejected after a close encounter between a binary and a single 
star. Therefore, mgys can be approximated to twice the aver¬ 
age mass of O stars in a cluster. We assume. 


^sys — 


2(to>17.5M0) = 2 


/ '^max /- / \ 1 

^ m^[m)dm 
17.5 Uw)dm 


777 

8.66^f^ 

rTTmax 


0.42 


(9) 


using the IME (Equation (1)). Thus nisys is a function of 
ttJmax SO that TUsys — m'sys(^max)- DuC tO the 777niax“.^ecl 
relation, irisys is thus a function of the cluster mass (the lower 
panel of Eigure 4). 

The number of ejections is proportional to the encounter 
rate, 

^‘’iVej^O OC r OC msys(TOmax)Mg°(;f. (10) 


The TTimax-^eci relation has only a minor effect on the num¬ 
ber of ejected O stars for the massive clusters, i.e., the depen¬ 
dence of nisys on Med becomes negligible in Equation (10). 
*^iVej^o is thus expected to be proportional to M^^f for the 
massive clusters. We fit (““Atgj o) as a power-law of cluster 
mass for clusters with 10^'^ < Med/M 0 < lO'^'^ (Eigure 3). 
These fits provide for each model sequence. 


(““iVej,o) CX 


' j^^0.63±0.066 
^0.78±0.034 

«^tfVeio.oei’ 
^0.62±0.061 

-'“eel , > 

» ^0.87±0.100 
-'“eel > 


MS30P, 

NMS30P, 

MS30P.SP, 

MS3UQ_SP, 

MS3S, 

MS80P, 

MSIOP. 


( 11 ) 


All model sequences but MS30P_SP and MSIOP exhibit 
slopes similar to each other within an error, in a range of 0.62- 
0.67. The power-law exponents we obtain from our A^-body 
models are slightly steeper than 0.5 in Equation (10) but with 
only a marginal difference. It can be argued that this is due to 
a broad range of binary parameters for which binaries could 
have engaged in the ejection of O stars in our models and the 
complexity of the few-body scattering process. The model 
sequence MS30P_SP shows a steeper slope, i.e., the more 
massive clusters in this model eject O stars more efficiently 
compared to other models. As the model sequence MSIOP 


begins with an extreme density the results of the model sig¬ 
nificantly deviate from the simplified theoretical expectation 
with simple approximations. 

Einally, the theoretically expected ejection fraction of O 
stars is 

th ^ 

*Vej.O = 77^^ CX F(m^ax) (12) 

7VO.IMF 


where F{m^sx) comes from the part of Equations (4) and 
(10) in which mmax plays a role. For massive clusters 
(Meci > 10^'® Mq), F(mniax) in Equation (12) becomes 
very weakly dependent on the cluster mass as in Equa¬ 
tion (10), and thus this term can be ignored. Therefore, **'/ej,o 
is expected to be proportional to in this cluster mass 

regime. Figure 5 shows the linear fits to the O star ejection 
fractions as a function of cluster mass for the massive clus¬ 
ters (10^'® < Meci/M© < lO"^-^) in a log-log scale. Within a 
cluster mass range lO^'^-lO"^'® Mq, the numerical data sug¬ 
gest the relation between O star ejection fraction and cluster 
mass as follows. 


r“/ej.0) CX 


7i^-0.40±0.064 
-'“ecL ’ 

,,;^-D.40±0.055 
-'“eel , ’ 

7i,e-D.25±0.029 
-'“eel , J 

;i,e-0.20±0.060 

-'“eel , J 

i^-0.43±0.055 
-'“eel , ’ 

a^-0.31±0.094 
-'“eel ’ 

I t-0.16±0.036 
-'“eel I 


MS30P, 

NMS30P, 

MS30P.SP, 

MS3UQ_SP, 

MS3S, 

MS80P, 

MSIOP. 


(13) 


The results of most models are reasonably consistent with the 
relation between the ejection fraction and cluster mass ex¬ 
pected from the binary-single star scattering encounter rate 
(Equation (12)) when applied to massive clusters in which the 
maximum stellar mass is less sensitive to cluster mass. 

For lower mass clusters, the number of O stars in a clus¬ 
ter and the ejection rate do not simply follow from the above 
equations due to the dependency of Wmax on the cluster mass 
and the highly collisional nature of the cluster cores which 
contain but a few O stars. 

The ejection fractions from individual clusters are plotted 
with open circles in Figure 2. For lower mass clusters the O 
star ejection fraction varies largely from one cluster to an¬ 
other, e.g., the values spread from 0 to 1 for clusters with 
Meci < 10^-® Mq. The spread gets smaller with increasing 
cluster mass as the stochastic effect diminishes with a larger 
number of O stars in more massive clusters. For massive clus¬ 
ters (Med > 10 "^ Mq) the individual ejection fractions mostly 
lie close to the average value. 

Only for a comparison, we additionally show results using 
10 pc as the ejection criterion in Figure 2 (without a constraint 
on the velocity of the stars) since a larger distance criterion is 
more appropriate from an observational point of view. The 
ejection fractions using 10 pc as a distance criterion show a 
similar shape and the peak at the same cluster mass, only the 
values are smaller (max. up to « 17 per cent for clusters with 
an initial rh = 0.3 pc) than the results using our main ejec¬ 
tion criteria. The difference becomes smaller or almost dis¬ 
appears for the massive clusters. Massive clusters may eject 
O stars most efficiently at earlier time of the evolution due 
to their short dynamical (i.e., crossing) time-scale, and soon 
their ejection efficiency would drop as the O star population is 
reduced in the core of the cluster. Also, in the massive clusters 
the velocity of the ejected systems is generally high (see Sec- 
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Figures. The average number of ejected 0-star systems, A^ej.o). 
a function of cluster mass. Lines are linear fits to ('^°“'Afej,o) with cluster 
mass between 10® ® and 10^ ® Mq (Equation (11)). 



2 3 4 5 


logio (J^ec./Mo) 

Figure 4. Top: Number of O stars in a cluster calculated using the IMF, 
-^OJMF (Equation (4), solid line). The dotted line indicates Nq oc M^d- 
Red filled circles present the average number of 0-star systems from A^-body 
calculations, in the MS30P model sequence. A red dashed line 

indicates {^°™Nq) oc obtained by a linear fit of logiQ{^°^No) - 

log ]^0 -^ecl to the MS30P model sequence. The other two model sequences, 
NMS30P and MS3S, have almost the same slopes as the MS30P model se¬ 
quence (Equation (5)). The grey solid line is the number of O stars derived 
from the IMF with the universal upper mass limit of 150 Mq. Bottom: As¬ 
sumed system mass of an O star binary, rrisys- The values are twice the av¬ 
erage O star mass in a cluster calculated from the IMF (Equation (1)). Here 
we assumed that an O star binary is an equal mass binary. Although the 
mass-ratio distribution of O star binaries is close to a uniform distribution for 
q > 0.1 in observations (Sana et al. 2012), our assumption is sufficient to 
show the dependence of the average system mass of O star binaries on cluster 
mass. 



Figure 5. The average ejection fraction of 0-star systems. Lines are linear 
fits of the average ejection fraction from the clusters with Med = 10^'^- 
10^'^ Mq (Equation (13)). The grey solid line shows the relation *^/ej,o oc 



Figure 6. Fitting function for (Equation (14) and Table 2). The 

symbols and line styles are the same as in Figure 3. 

tion 5.1) so they quickly move further than 10 pc away from 
the cluster centre. Low-mass clusters have longer dynamical 
times-scales and they eject O stars moving relatively slowly 
so that they need more time to reach a distance of 10 pc. 

4. FIELD O STARS FROM DYNAMICAL EJECTION PROCESSES 

Our A^-body results show that logio(^°™-^ej,o) steeply in¬ 
creases with logio -^eci ^ lower cluster-mass range (Med < 
10^'^ Mq) while it can be approximated to a linear func¬ 
tion of logio -^eci at a higher cluster-mass range. In the 
previous section, we use data only from the clusters with 
10^'^ Mq < Med < 10^'^ Mq to fit a linear relation between 
logio(‘^°”^-/Vej,o) and logio -^eci (Figure 3). In order to fit the 
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Table 2 

Parameters by fitting Equation (14) for all models. 


Model 

a 



b 



n 


MS30P 

-0.56 

± 

0.25 

0.53 

± 

0.07 

0.29 

± 

0.04 

NMS30P 

-0.43 

± 

0.25 

0.49 

± 

0.07 

0.42 

± 

0.05 

MS30P.SP 

-0.91 

± 

0.16 

0.64 

± 

0.04 

0.26 

± 

0.04 

MS3UQ.SP 

-1.04 

± 

0.25 

0.67 

± 

0.07 

0.22 

± 

0,04 

MS3S 

0.67; 

tO.32 

0.14 

± 

0.08 

0.90 

± 

0.10 

MS80P 

-0.60 

± 

0.01 

0 

CO 

-,3 

±1 

0.003 

0.51 

± 

0.03 

MS I OP 

-1.17 

± 

0.16 

0.80 

± 

0.05 

0.16 

± 

0.03 


number of ejected O stars for the entire cluster mass range 
in which is > 0, we use the following functional 

form, 


logio y 


a + 6 logio ^ ~ 


1 

(logio a: - logio 


(14) 


where x > xq, y = (“™A^ej.o) and x = M^ci/Mq. In this 
function, logio V converges to — oo for x ^ xq and approx¬ 
imates a linear function, a -f & logio x, x ^ xq. This 
function can thus be fitted to our results. The value xq in¬ 
dicates the cluster mass where the y value, i.e., (“™7Vej_o)^ 
Ri 0. For all models but MS3S and MS80P, we choose 
Xq = 240 Mq which is the maximum cluster mass at which 
clusters do not host single O stars initially, as given by the 
Wmax-lUeci relation. While for the other two models we 
choose xo = 10^ ® Mq at which mass a cluster contains O 
stars but (““A^ej,o) is 0 in the N-body calculations. By per¬ 
forming non-linear least squares fitting of Equation (14) to our 
results, we obtain the parameters listed in Table 2 and plot the 
results in the upper panel of Figure 6. 

The mass distribution of young star clusters can be de¬ 
scribed with a simple power-law. 


^ecl(Mecl) = fceclM-f, (15) 

where the normalisation constant feed is calculated from the 
total mass of all clusters which form within a certain time. 
Assuming a global Milky Way star formation rate (SFR) of 
3 Mq yr“^, a minimum cluster mass Med,min = 10 Mq, 
and P = 2 for the Milky Way, the maximum cluster mass, 
Med,max, that Can form in the Galaxy becomes Ri 1.9 x 
10® Mq (Weidner et al. 2004; Marks & Kroupa 2011, and 
references therein). The total mass formed in stars over time 
St, Mscs, is 


Mscs = SFR x6t = 


pMecl , max 
> Afecl ,min 


Med Cecl(Med) rfMed 

(16) 


where 5t is the cluster-population formation time-scale, about 
10 Myr (Weidner et al. 2004). We adopt here the Med,max- 
SFR relation derived by Weidner et al. (2004), i.e., the mass 
of the most massive forming cluster is dependent on the phys¬ 
ical properties of the star-formation environment provided by 
a self-regulated galaxy. It has been found that the luminosity 
of the brightest cluster increases with the galaxy-wide SFR 
(Farsen 2002, and references therein) or with the total number 
of clusters (Whitmore 2003). The relation is treated to be the 
result of the size of the sample in several studies (Whitmore 
2003; see also section 2.4 of Portegies Zwart et al. 2010, and 
references therein). Recently, Pflamm-Altenburg et al. (2013) 
ruled out that the relation results from the pure size-of-sample 
effect by studying a radial dependancy of maximum star clus¬ 
ter masses in M33. Randriamanakoto et al. (2013) confirm 


X 


10" 



10 " 


10“* 


10 " 


(M®) 



Figure 7. Top: Total number of the ejected O star systems per cluster mass 
in the Galaxy (Equation (17)). For example, according to the most realistic 
model sequence (MS3UQ_SP), in total about 4 ejected O star systems per 10 
Myr are contributed from all clusters weighing about 400 Mq. The small 
numbers of ejected O star systems contributed from massive clusters come 
about because the number (or probability) of such massive clusters to form is 
small according to the very young or embedded cluster mass function (Equa¬ 
tion (15)). Bottom: Cumulative number of the ejected O star systems as a 
function of cluster mass. Right y-axis shows the ejected O star fraction of 
the total O stars formed within 10 Myr deduced by using Milky Way type 
parameters. Numbers next to the model sequence names are the expected 
total number of ejected O star systems which are the sum of Equation (17) 
over all clusters for each model sequence. For example, the MSlOP model 
sequence would imply that all clusters formed within 10 Myr eject in total 
almost 2.5 x 10^ O star systems assuming the SFR = 3 Mq yr~^ such 
that Med max = 1-9 10^ Mq. As another example, according to the most 

realistic sequence of models (MS3UQ_SP), every 10 Myr, 10818 O star sys¬ 
tems are ejected from all clusters formed in this time interval, clusters with 
a mass of almost 400 Mq have the largest contribution to this population 
(upper panel) and clusters with masses up to 10^ Mq contribute almost 33 
per cent to this population (lower panel, right axis), 50 per cent of all ejected 
O stars originate from Med < 3 x 10^ Mq clusters. 


this by studying the Med,max versus SFR relation in their in¬ 
dependent observational survey. The authors noted that their 
result can be explained with purely random sampling but only 
if the luminosity function of super star clusters at the bright 
end is very steep (a slope > 2.5), steeper than usually ob¬ 
served 2). 

By combining the fitting function for the number of ejected 
O stars (Equation (14)) with the cluster mass function (Equa¬ 
tion (15)), the number of O stars being ejected to the field con¬ 
tributed by clusters in the mass range -Med lo Med + d^eci 
is 

d‘°*"*iVej.o(Mecl) = (™“(Vej,o)(Mecl) X CecKMed) dMed, 

(17) 

and is shown in the upper panel of Figure 7. As the number of 
clusters decreases more steeply (^eci(Med) oc M“j^) than the 
number of ejected O star systems from a cluster increases with 
cluster mass ((“™Vej,o)(Med) oc Mj^j with 0.6 < 7 < 0.9, 
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Section 3), the number of O stars contributed to the field by 
dynamical ejections is a function decreasing with increasing 
cluster mass. The most prominent contributors to the popula¬ 
tion of ejected O stars are f«400 Mq clusters. 

In the lower panel of Figure 7, cumulative numbers of 
ejected O-type systems as a function of cluster mass are 
shown. The fraction (right Y-axis) is calculated by dividing 
the cumulative number by the total number of O stars formed 
during the cluster-population forming time scale of 10 Myr, 
which can be deduced from Equations (4) and (15), 


total 


Atecl ,max 
^ec\ ,min 


-/Vo,IMF(-M^ecl) ^ecl(Afecl) dMgcl- 


( 18 ) 

The parameters which we adopt here for the Milky Way 
provide totaijy^ 66749 being formed every 10 Myr. 
Note total(.jjg number of individual O stars 

which have formed within the formation time scale, while 
d*°taiWej Q(Meci) is the number of ejected O star systems. It 
is difficult to derive the total number of O star systems from 
the models because a binary fraction evolves due to dynami¬ 
cal disruptions/captures and mergers between stars, the final 
binary fraction varies with cluster mass and initial conditions, 
and some O stars have an O star secondary while some do not. 
The fraction shown in the lower panel of Figure 7 can, thus, be 
interpreted as being a minimum value since the total number 
of O star systems would be smaller than which counts 

all individual O stars from the IMF and cluster mass function. 
Our models thus yield as a lower boundary that binary-rich 
clusters with an initial rh = 0.3 pc can eject Ril3-16 per cent 
of the total O star systems and even Ri 38 per cent can be dy¬ 
namically ejected if clusters form more compact, e.g., with an 
initial rh = 0.1 pc (Figure 7). About 50 per cent of ejected O 
stars come from clusters with Med ^ 3000-5000 Mq. 

The fraction of O stars found in clusters and/or OB- 
associations is 70-76 per cent (de Wit et al. 2005; Eldridge 
et al. 2011, depending on how to count a binary) yielding the 
fraction of the field O stars as being 24-30 per cent. Our 
model is thus reasonably consistent with these data. Addi¬ 
tional three processes that are not part of the present model 
but which increase the apparent number of field O stars are : 
(i) Some of the observed field stars may originate from clus¬ 
ters which contained only a few or one O stars and which dis¬ 
solved at an early age by rapid gas removal (Gvaramadze et al. 
2012). That is such O stars “lost” their birth cluster through 
disruptive gas removal (see also Kroupa & Body 2002; Weid- 
ner et al. 2007). (ii) Small-iV clusters for which the two-body 
relaxation time is comparable to the crossing time also dis¬ 
solve within a few initial crossing times (Moeckel et al. 2012), 
which is aided by stellar evolution within a few tens of Myr. 
(iii) Some O stars may be ejected when they are companions 
to primaries that explode as supernovae at an age older than 
3 Myr (Eldridge et al. 2011). 

In observations about 6 per cent of O stars are found be¬ 
ing runaways which have a velocity higher 30 kms“^ (El¬ 
dridge et al. 2011). The runaway fraction among the dynam¬ 
ically ejected O stars generally increases with cluster mass 
(upper panel of Figure 8), e.g., 0 (Med = 10^'^ Mq) to 43 
(Med = 10"^ ® Mq) per cent in the case of the MS3UQ_SP 
model and 81 percent for the 10^ Mq model in Banerjee et al. 
(2012). We deduce the number of the runaways first by fitting 
a function for the runaway fraction as a function of cluster 
mass with the same functional form of the number of ejected 



■^ecl (M©) 


Figure 8. Top: Runaway fraction among the ejected O stars from two model 
sequences, MS30P_SP and MS3UQ_SP. The data are fitted with the same 
functional form of Equation (14) (curved lines). Error bars indicate Poisson 
errors. Bottom: Expected cumulative number of runaways as a function of 
cluster mass. The numbers are derived by multiplying the fitted lines in the 
upper panel to the number of ejected O star systems with cluster mass, Equa¬ 
tion (17). More than 70 per cent of O star runaways are expected to originate 
in the clusters more massive than 10^ AT©. 

O stars (Equation (14), but y being the runaway fraction) 
and then by multiplying the fitting function by the number of 
ejected O stars (Equation (17), upper panel of Figure 7). By 
dividing the total number of the runaways by the total num¬ 
ber of the ejected stars (lower panel of Figure 7), we obtain 
the runaway fraction among the ejected O stars. Our most 
realistic model, MS3UQ_SP, predicts about 30 per cent of the 
ejected O stars, i.e., about 5 per cent of all O star systems to be 
runaways. Eldridge et al. (2011) studied runaways assuming 
the binary supernova scenario, in which the companion star 
obtains a high velocity when one of the binary components 
undergoes a supernova explosion. They predicted about 0.5 
per cent ^ of O stars being runaways when only this process is 
accounted for. The value combining ours and Eldridge et al. 
(2011) is thus well consistent with the observed one. 

Figure 7 demonstrates that the presence or absence of ini¬ 
tial mass-segregation leads to differences in the numbers of 
ejected O stars only for lower mass clusters (Mgd < 10^ Mq) 
with the constant half-mass radius we assume in this study. If 
larger rh are assumed the cluster mass at which the differ¬ 
ence vanishes would move up to a larger cluster mass since 
the mass-segregation time-scale gets longer. 

5. PROPERTIES OF EJECTED O STAR SYSTEMS 

In this section, we study how some properties, such as ve¬ 
locities, distances from the cluster centre, and system masses 

^ Eldridge et al. (2011) provides two values, 0.5 (mentioned in the text) 
and 2.2 per cent, for the fraction of 0-star runaways produced by a supernova 
explosion of the primary star in a binary system. The former is obtained 
using 0 stars with velocity larger than 30 km s“ ^, which is also our runaway 
criterion, the latter is derived with a lower velocity limit of 5 kms“^ for 
runaways. The latter value is applicable to process (iii) in the text. 
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of the ejected O-star systems vary with cluster mass. We 
remind the reader that our analysis is restricted to systems 
younger than 3 Myr. 

5.1. Kinematics of the ejected O star systems 

Cumulative distributions of velocity, distance from the clus¬ 
ter centre, and system-mass of all ejected O-star systems for 
each cluster mass and for each set of cluster models are shown 
in Figure 9. Over all, the systems ejected from more massive 
clusters have velocity distributions skewed towards higher ve¬ 
locity. This can be understood in two ways. Firstly, the veloc¬ 
ity dispersion of stars in a cluster increases with cluster mass 
under the assumption of a constant cluster size (in general, ve¬ 
locity dispersion is a function of both size and mass of a clus¬ 
ter), i.e., on average stars in more massive clusters have higher 
velocities. Secondly, a star needs a velocity greater than the 
escape velocity which also increases with cluster mass, in or¬ 
der to leave its parental cluster. For the MS80P model se¬ 
quence (with Th = 0.8 pc), the trend is not as pronounced 
as in other models since the ejection fractions are small and 
so are the total number of ejected O-star systems, e.g., only 
in total no stars are ejected out of 4 runs in the Mq 
cluster models. Also in the case of the lO^’® Mq clusters in 
the NMS30P models (rh = 0.3 pc) only three O star systems 
are ejected in 100 runs. 

The highest velocity of an ejected system generally in¬ 
creases with cluster mass, although this is not always the 
case, e.g, in the MS30P model sequence, the system mov¬ 
ing the fastest (a single star of mass 46 Mq with a velocity 
of « 200 km s“^) is ejected from a 10^’® Mq cluster. The 
velocity distribution is dependent on the cluster model, es¬ 
pecially on the initial half-mass radius of the cluster. The 
MSI OP models have velocity distributions biased towards 
higher velocities compared to the larger-sized cluster mod¬ 
els (Figure 9). The highest velocity of an ejected system (a 
single star of mass 34 Mq from an Med = lO"*’® Mq cluster) 
is « 300 km in the MSI OP model sequence. 

As the distance distribution of the ejected systems is di¬ 
rectly related to the velocity distribution, their shapes are sim¬ 
ilar. The slight difference between them results from the time 
of ejections occurring at systematically different times in dif¬ 
ferent models. Owing to the systematic shift to larger ejection 
velocities, the ejected systems travel on average further away 
from the clusters in the more massive clusters. The distances 
of the ejected systems from Med > 10'^ Mq clusters ranges 
from a few pc to a few 100 pc. The most extreme system (a 
single star with a mass of 34 Mq) ejected from an MSlOP 
model traveled up to « 750 pc away from its birth cluster by 
3 Myr. 

Cumulative mass distributions of the ejected O stars are 
plotted in the right column of Figure 9. In the case of a binary 
system, masses of both components are counted separately. 
As binary fractions of the ejected systems are low for most 
of the models, the cumulative distributions of system masses 
would be similar to those of individual stellar masses. A few 
models suggest that the more massive clusters may have a 
higher fraction of more-massive ejected stars. However the 
trend is not strongly evident. For example, in the case of the 
MSlOP models, the mass distributions of the ejected O stars 
look similar to each other for different cluster masses. But 
only clusters with Mgd > 10^'^ Mq eject O stars heavier than 
100 Mq as the lower mass clusters do not have such a massive 
star initially, e.g., for Mgd = 1000 Mq, rrimax ~ 44 Mq. 
About 15 per cent of the ejected O stars from clusters with 


-^eci ^ 10^’® Mq in the MS30P_SP models are more mas¬ 
sive than 100 Mq, while this is the case for < 5 per cent from 
the other models. This may result from a higher number of 
mergers in the MS30P_SP models due to the high fraction of 
close, eccentric binaries in comparison to the other models. 

Thin grey lines in the right column of Figure 9 indicate 
cumulative mass distributions of O stars from the canonical 
IMF with the lower- and the upper-mass limit of, respec¬ 
tively, 17.5 Mq and the mass of the most massive star among 
the ejected O star for each cluster-mass set. It follows that 
most models show that ejected stars have a shallower slope of 
the mass function than the canonical value (index of 2.3, see 
Equation (1)) as also shown in Fujii & Portegies Zwart (2011) 
resembling the mass function in the cluster core. 

The normalised cumulative number of the ejected O star- 
systems as a function of time for a few models are shown 
in Figure 10. In the case of the MS30P_SP model (Mgd = 
10^'^ Mq), almost 70 per cent of the ejected stars travel be¬ 
yond 3 times the half-mass radius before 1.5 Myr. Initially 
not mass-segregated cluster models (e.g., the 10^'® Mq clus¬ 
ters of the NMS30P model sequence in Figure 10) show a 
delay of the ejections as the model requires time for massive 
stars to migrate to the cluster centre to interact with other mas¬ 
sive stars. But about half of the ejected stars already reach 
a distance of 3 times the half-mass radius from the cluster 
centre before 1.5 Myr of evolution. Also the ejection rate is 
higher during the earlier times of evolution. Note that the 
time at which a star is classified as being ejected is differ¬ 
ent to when the ejection has actually happened. The time 
presented in the figure is comprised of the traveling time the 
stars need to reach to the ejection criterion and thus the ejec¬ 
tions actually happened at an earlier time than shown in the 
figure. This can naturally explain the steeper increase of the 
number of ejected stars at a later time for the Banerjee et al. 
(2012) model in which the initial half-mass radius of the clus¬ 
ters is larger (rh = 0.8 pc) than in the above two models 
(rh = 0.3 pc). Thus an ejected star requires more time to 
reach the ejection criterion in their model (see Section 3 for 
the criteria). A rough estimate of the time when the ejections 
occur peaks at around 1 Myr, though it can vary with different 
initial conditions (Oh et al. in preparation). 

5.2. The binary fraction among the ejected O star systems 

Not only single O stars but binaries containing O star com¬ 
ponent/s are ejected from star clusters via dynamical interac¬ 
tions. The binary fraction of ejected O stars can be substan¬ 
tial, up to « 78 per cent (Table Al). Even initially single-star 
clusters eject O star binaries dynamically (Table Al). In Eig- 
ure 11 is plotted the average binary fraction among ejected 
O-star systems from the binary-rich cluster models as a func¬ 
tion of cluster mass. A smaller fraction of binaries is ejected 
from the more massive clusters because the kick required to 
remove a system from the cluster is larger and more binaries 
are disrupted through interactions with other stars in the more 
massive clusters. 

The binary fraction among ejected O stars can be signifi¬ 
cantly altered by assuming a different initial period distribu¬ 
tion. An initially higher fraction of close binaries leads to 
a higher binary fraction of ejected O stars. Thus the more 
realistic MS30P_SP and MS3UQ_SP models produce higher 
binary fractions of ejected O star systems compared to the 
other models (Eigure 11). The observed spectroscopic binary 
fraction among the Galactic field O stars is 43 - 50 per cent 
(Gies 1987; Mason et al. 1998, 2009; Chini et al. 2012). This 
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Figure 9. Normalised cumulative distributions of velocity {left) and of distance from the cluster centre {centre) of ejected O star systems, and of mass of the 
individual ejected O stars {right) from our calculations. See the text for the grey lines in the right most column. 
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Figure 9. (Continued) 
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Figure 10. Cumulative O star ejection as a function of time for 10® ’® Mq 
clusters in models MS30P_SP and NMS30P, and the Banerjee et al. (2012) 
model. For example, almost 45 per cent of all ejected 0 star systems are 
ejected by 1 Myr in the MS30P_SP (10® ® Mq) model. 
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Figure 11. Top: Averaged binary fraction among the ejected 0-star sys¬ 
tems, (/b,ej,o)’ froni initially binary-rich cluster models. A (grey) cross at 
Meci = 10^ Mq is the Banerjee et al. (2012) model which adopted a uni¬ 
form initial period distribution for O star binaries. Error bars indicate stan¬ 
dard deviations of (/b,ej,o)- The observed spectroscopic binary fractions in 
the Galactic field O stars are shown as a grey area (Gies 1987; Mason et al. 
1998, 2009; Chini et al. 2012). Bottom: Binary fraction among the runaway 
O stars from initially binaiy-rich cluster models. The cluster models that do 
not eject any runaway O star, e.g., 10^'^ Mq cluster models from all se¬ 
quences but the MSlOP sequence, are not included in the figure. En'or bars 
indicate Poisson errors. The grey area indicates a range of the observed spec¬ 
troscopic binaiy fractions of runaway O stars in Gies (1987) and Mason et al. 
(1998, 2009). The grey horizontal line is the observed value in Chini et al. 
( 2012 ). 
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may be reproduced readily through dynamical ejections of 
O stars with more-realistic (than the Kroupa 1995b) period 
distribution function for O star binaries, such as constrained 
recently by Sana et al. (2012) (Equation (3)). This comes 
about because, as can be seen for the most realistic sequence 
(MS3UQ_SP), the dominant contributors to the population of 
ejected O star systems are < 10"^ Mq clusters (lower panel 
of Figure 7), and these ejected O star systems have a binary 
fraction > 30 per cent (upper panel of Figure 11). It should be 
noted that the given observed binary fractions are lower limits 
of the binary fraction of the held O stars since only spectro¬ 
scopic binaries are counted. 

In the case of runaways, the binary fractions from our most 
realistic sequences, MS30P_SP and MS3UQ_SP, are well 
consistent with the observed spectroscopic binary fraction of 
the runaway O stars in Gies (1987) and Mason et al. (1998, 
2009), 19-29 per cent (grey area in the lower panel of Fig¬ 
ure 11), but too small compared to the value in Chini et al. 
(2012),^ 69 per cent (grey horizontal line in the same hgure). 
This may be due to the difference in assigning runaways be¬ 
tween the Chini et al. (2012) work and the other studies. In 
Chini et al. (2012), stars are classihed as runaways when they 
can be traced back to any cluster or associations in the study 
of Schilbach & Roser (2008). In other words, the runaways in 
their sample may be simply ejected O stars from star clusters 
which have masses < 10^ Mq since very massive clusters are 
exceedingly rare. The other studies assign stars to runaways 
with a large distance from the Galactic plane or with a high 
peculiar (radial) velocity. 

In most of our cluster models (with the exception of the 
MS30P_SP and MS3UQ_SP models) the binary fraction of 
O stars, whether they are inside a cluster or not, at 3 Myr 
is smaller than the observed O star binary fraction (« 70 
per cent, Chini et al. 2012; Sana et al. 2012, and references 
therein). This is the case because dynamical interactions be¬ 
tween binaries and other cluster members disrupt binaries, 
preferentially those with a wide orbit (Kroupa 1995a; Parker 
et al. 2009; Marks et al. 2011) and because the Kroupa pe¬ 
riod distribution (which was constrained originally only for 
late-type stars) extrapolated in some of the model sequences 
to the O star regime has a larger fraction of long period bina¬ 
ries compared to the observed distribution of O star binaries 
(Figure 1). 


6 . DISCUSSION 

While the mass of clusters which eject O stars most effi¬ 
ciently is the same for all models (Med ~ Mq), it 

has been shown that the dynamical ejection of O stars can 
be highly sensitive to the initial conditions of star clusters. In 
this section, we further discuss the effects of our choice of ini¬ 
tial parameters for massive stars and for the clusters, and the 
available observational constraints. 

The ejection fraction is strongly dependent on the initial 
cluster radius. As seen in Figure 2 the shapes of the curves 
look similar but the values of the ejection fraction are larger 
for initially smaller sized clusters. Thus it is important to 

® We note here that the Chini et al. (2012) sample is not biased. They 
selected stars from the Galactic 0 Star Catalogue v.2.0 (Sota et al. 2008) of 
which 248 stars could be reached from their observatory. They additionally 
employed archival data from the European Southern Observatory (ESO) for 
those stars already contained in their survey. Thus the ESO spectra did not 
add any new star to their sample but increased only the number of spectra 
and, as a consequence, the time basis for individual stars (R. Chini, private 
communication). 
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check if our models are comparable to observed young star 
clusters. The half-mass densities of our iV-body models at 
3 Myr range from 175 to « 10® Mqpc~^, from « 160 to 
« 3.4 X 10^ M 0 pc“®, and from « 50 to « 4000 Mqpc~^ 
for initial = 0.1, 0.3, and 0.8 pc clusters, respectively. 
These values are well within the densities of observed young 
(< 5 Myr) massive (Med > 6000 Mq) star clusters, clas¬ 
sified as star burst clusters in Pfalzner (2009), which range 
from 1500 to 4 x 10® Mqpc~^ (Figer 2008). However, cau¬ 
tion is due in directly comparing the model cluster densities to 
the observed values since the observed densities are closer to 
the central (i.e., highest) densities rather than to the half-mass 
density. Concerning the initial radii of the embedded clusters, 
these are dynamically sub-pc (Testi et al. 1999; Lada & Lada 
2003; Kroupa 2005; Marks & Kroupa 2011; Smith et al. 2005; 
Tapia et al. 2009, 2011, 2014). 

The initial binary population plays a major role in the dy¬ 
namical ejections of massive stars (see the differences be¬ 
tween the MS30P and MS3S models). Thus it is important 
to adopt realistic initial binary parameters to understand mas¬ 
sive star ejections in reality. The pairing method for mas¬ 
sive binaries adopted in most of our models results in extreme 
mass-ratios, e.g., most of them are larger than 0.9. This was 
motivated by the high fraction of twins in observed massive 
binaries (Pinsonneault & Stanek 2006; Kobulnicky & Fryer 
2007, but see also Lucy 2006), while recent observational 
studies favour a rather flat mass-ratio distribution for O-type 
binaries (Kiminki & Kobulnicky 2012; Sana et al. 2012) still 
excluding random pairing from the IMF for this spectral type. 
Adopting a more realistic initial mass-ratio distribution for O- 
type binaries may reduce the efficiency of the ejection, though 
the shape of the /ej^o(-^eci) ejection fraction as a function 
of parent birth-cluster mass, is expected to be kept. In this 
contribution we include a realistic model for the observed pe¬ 
riod distribution of O star binaries (Equation (3)). The model 
with the observed period distribution, which reproduces the 
observed abundance of short period binaries, exhibits a sub¬ 
stantial increase of the ejection fraction, by « 30 per cent, 
from massive clusters (Med > 10'^ Mq). Whether the ob¬ 
served distribution is a real initial distribution is still unclear 
as the observed distribution has been compiled from several 
clusters which may have already evolved dynamically despite 
their young ages (1^ Myr). The authors suggested that dy¬ 
namical evolution would be negligible for such massive close 
binaries (log]^Q(P/days) < 3.5). However, the population of 
massive binaries can evolve through the dynamical interac¬ 
tions with other massive stars/binaries in the clusters (Oh et 
al. in preparation). Moreover, short-period massive binaries 
are vulnerable for merging either due to stellar evolution or 
the perturbation-induced collision of the two binary-star com¬ 
panions. Thus the effect of dynamical evolution needs to be 
taken in to account iteratively to find the true initial period 
distribution of massive binaries. But this is beyond the aim 
of this study. For late-type binaries this has already been per¬ 
formed successfully by Kroupa (1995a,b, see also Marks et al. 
2011; Marks & Kroupa 2012). 

Unlike the other initial conditions mentioned above, initial 
mass segregation does not result in a difference of the O star 
ejection fraction for clusters with Med > 10 ® Mq and ini¬ 
tially rh = 0.3 pc due to efficient dynamical mass segrega¬ 
tion. But differences do occur for lowest mass clusters: ini¬ 
tially mass-segregated clusters more efficiently eject O stars 
than non-segregated clusters, as dynamical mass segregation 
is inefficient to gather massive stars in the core of such clus¬ 


ters. For clusters with larger initial rh values, as the dynamical 
mass segregation time scale becomes longer, the initial mass 
segregation can play a significant role. In observations, post 
gas-expulsion young massive clusters (age < 5 Myr) are gen¬ 
erally compact (effective radii < a few pc, Portegies Zwart 
et al. 2010; Kroupa 2005; Tapia et al. 2003). 

Throughout this study the mmax-ATeci relation is adopted. 
It follows from an intuitive theoretical concept that the mass 
of a star forming molecular cloud core is finite and that 
star formation is a self-regulated process (Adams & Fatuzzo 
1996), and it follows from the data on young stellar popula¬ 
tions (Weidner et al. 2010, 2013a; Kirk & Myers 2011). The 
existence of a physical Wniax-Afeci relation is thus taken to 
be established (Weidner et al. 2014). Although the theory 
has been successful in describing the stellar contents of whole 
galaxies (e.g., Weidner et al. 2013b), it has been debated if the 
relation exists (see section 2.1 in Weidner et al. 2013a, for de¬ 
tails) or if star formation is a purely stochastic process (e.g., 
Krumholz 2015). This is a fundamental issue for understand¬ 
ing galaxy-wide IMFs and how massive stars form, and we 
here thus stress that the data do not support the hypothesis that 
star formation is stochastic (see section 7 of Kroupa 2015). 
For example, Hsu et al. (2012, 2013) presented that LI641, 
the low-density star-forming region (containing as many as 
« 1600 stars) of the Orion A cloud, is deficient of massive 
(O and early B) stars compared to the canonical IMF with 
3^(7 significance and that with a probability of only 3 per 
cent the southern region of L1461 and the ONC can be drawn 
from the same population supporting that the high-mass end 
of the IMF is dependent on environmental density. In sec¬ 
tion 4.2 of Krumholz (2015) three references (Calzetti et al. 
2010; Fumagalli et al. 2011; Andrews et al. 2013) are cited as 
containing evidence against a physical mmax-Med relation. 
Weidner et al. (2014) critically discuss these works showing 
Andrews et al. (2013) data to be entirely consistent with the 
physical mmax-Med relation. As already discussed in Sec¬ 
tion 3, whether the mmax-Med relation is adopted or not does 
not affect the here obtained results though. 

We do not include a gas potential, which would subse¬ 
quently be removed rapidly, in the calculation (cf., Kroupa 
et al. 2001; Banerjee & Kroupa 2014). The main outcomes of 
gas expulsion are significant stellar-loss from and expansion 
of clusters and their possible subsequent dissolution. There¬ 
fore gas expulsion can weaken the ejection efficiency by low¬ 
ering the (central) density of a cluster if that occurs at a similar 
time when a majority of the ejection processes take place. In 
the future, further studies including a gas potential are needed 
to understand the effect of gas expulsion on the ejection of 
massive stars. Inclusion of a gas potential will be possible 
once the statistical quantification of the gas expulsion process 
is in place and computable (cf., Banerjee & Kroupa 2014). It 
has been argued recently that gas expulsion, if the process oc¬ 
curs at all, may have no effect on star cluster evolution (e.g., 
see Longmore et al. 2014, and references therein). It may, 
however, be difficult to detect the expanding signature result¬ 
ing from residual gas expulsion since revirialization after gas 
expulsion can be rapid for massive clusters, Med >10® Mq 
(Portegies Zwart et al. 2010; Banerjee (& Kroupa 2013). 

To constrain which initial conditions are responsible for 
the observed field O star populations, one needs the observed 
quantities to compare with the results from models. However, 
it is very difficult to find which clusters individual field O stars 
originate from. Some ejected O stars can travel further than a 
few hundred pc away from their birth cluster within a few Myr 
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thus, without knowledge of their full-3D kinematics and the 
Galactic potential, it is almost impossible to find their birth 
cluster. Gaia may help to find the birth clusters of the field 
O stars. However, the two-step ejection process of Pflamm- 
Altenburg & Kroupa (2010) complicates this. Their calcula¬ 
tion suggests that 1^.6 per cent of O stars would appear to 
have formed in isolation. They used the observed values of 
the O star runaway fraction (10-25 per cent, Gies & Bolton 
1986; 46 per cent. Stone 1991), and the binary fraction among 
the runaways (10 per cent, Gies & Bolton 1986). Since most, 
if not all, ejected O-star-O-star binaries may go through the 
two-step ejection process and then cannot be traced back to 
their birth cluster, the fraction of O stars which appear to have 
formed in isolation can be higher if a binary fraction of the 
dynamically ejected stars is substantial. A quantitative study 
is required to determine the fraction of O-stars experiencing 
the two-step ejection process. 

7. SUMMARY 

We study how the ejection of O-stars varies with cluster 
mass (or density as we assume a constant size for different 
cluster masses) under diverse initial conditions. We find that 
clusters with a mass of lO^’® Mq most efficiently shoot out 
their O stars to the field compared to other ones with other 
masses, for all models, i.e., this result obtains independently 
of radius or the presence of initial mass segregation or bi¬ 
naries. Our results suggest that moderately massive clusters 
(Med = 10^-^ Mq), which have formed about 10 O star sys¬ 
tems, most efficiently eject their O stars through energetic en¬ 
counters between massive stars in the cluster core. Up to on 
average « 25 per cent of the initial O star content is ejected 
from the clusters with this cluster mass in rh = 0.3 pc models 
(52 per cent in rh = 0.1 pc clusters. Figure 2). However, the 
spread is large and ejections of all O stars leading to remnant 
young clusters void of O stars occur in about 1 per cent of all 
such clusters (Figure 2). Fujii & Portegies Zwart (2011) sug¬ 
gested that one bully binary, dynamically formed, is mainly 
responsible for flinging out the stars from initially single-star 
clusters. Here we show that the number of ejections is sig¬ 
nificantly higher when the highly significant initial massive 
binary population is accounted for. The ejections, thus, oc¬ 
cur via close encounters involving several binaries in reality 
since the massive stars are observed to have high multiplic¬ 
ity fractions in young star clusters (Chini et al. 2012; Sana 
et al. 2012, and references therein). Furthermore, more realis¬ 
tic clusters, which have O-star binaries with periods following 
the Sana et al. (2012) distribution (MS30P_SP), lead to even 
larger ejection fractions, especially for massive clusters. The 
mass-position of the peak of /ej,o = /ej,o(-Meci), however, 
does not change (Med,peak « 10^'^ Mq). 

The decrease of the O star ejection fraction with increasing 
cluster mass for massive (Med > 10^'^ Mq) clusters can be 
represented by a simple power-law of cluster mass with an 
exponent between —0.16 and —0.43. Clusters with an initial 
half-mass radius rh = 0.3 pc and with relatively simple initial 
conditions (Kroupa initial period distribution for O stars or no 
initial binaries) have an exponent of « —0.4 being in good 
agreement with —0.5 which is derived from the number of O 
stars from the IMF and the binary-single star scattering rate 
(Equation (12)). Models with a high fraction of close initial 
binaries or smaller rh deviate from the expected theoretical 
value by having larger ejection fractions. 

We show that not only the ejection fraction but also the 
properties of the ejected O star systems depend on the initial 


cluster mass. The velocity distribution of the ejected O star 
systems shifts to higher velocities for more massive clusters, 
and so does their distance distribution from the cluster. But 
the distribution of ejected system masses varies weakly with 
cluster mass. The binary fraction of the ejected systems de¬ 
creases with increasing cluster mass as does the binary frac¬ 
tion of the systems remaining in the cluster because binary 
disruption via interactions with other members of the cluster 
is more efficient for a denser cluster, i.e., more massive clus¬ 
ter. The binary fraction among the ejected O star systems is 
substantial, especially in the more realistic models in which 
the period distribution of the observed O star binaries is im¬ 
plemented (MS30P_SP and MS3UQ.SP). 

Considering that the cluster mass function is approximately 
a power-law with a near Salpeter index, the main fraction of 
field O stars dynamically ejected from clusters originate from 
low/intermediate mass clusters (Mgd < a few 1000 Mq, see 
Section 5). The dynamical ejection process populates « 16 to 
38 per cent of all O stars, formed in clusters, to the field. Com¬ 
bining our results on the dynamical ejection fraction and the 
binary-supernova scenario computed by Eldridge et al. (2011) 
it follows that the observed fractions of the Galactic field- and 
runaway O stars are well accounted for theoretically. It should 
be noted that very massive (> 100 Mq) runaways can only 
originate from massive star-burst clusters (Mgd > 10^ Mq) 
because only massive clusters can harbour a sufficient number 
of very massive stars to eject some of them. 

We are grateful to Sverre Aarseth for making his NBODY6 
code freely available and for continuing its improvements. We 
also thank Rolf Chini for comments on the observational data 
in his paper. 
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APPENDIX 

A TABLE FOR THE RESULTS OF INDIVIDUAL MODELS 

Here we append the table containing the results and properties of each model cluster (Table Al). 
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Table A1 

Results of all models at 3 Myr. 


Med (Mq) 

^^max (M^0) 

(j-h) (pc) 

(comlVo) 

r“(Vej,0> 

/com f ^ \ 

\ /ej,0/ 

{/b,o) 

(/b,ej,o) 

A^run 

MS30P 

102.5 

21.2 

0.60 

1.11 

0.09 

0.080 

0.84 

0.78 

100 

10 ^ 

43.9 

0.72 

2.96 

0.61 

0.172 

0.66 

0.33 

100 

103.5 

79.2 

0.70 

9.57 

2.38 

0.252 

0.45 

0.17 

100 

10 ^ 

114.7 

0.50 

30.10 

4.40 

0.149 

0.23 

0.13 

10 

104.5 

136.2 

0.50 

105.00 

11.25 

0.105 

0.15 

0.05 

4 

NMS30P 

102.5 

21.2 

0.56 

1.14 

0.03 

0.020 

0.87 

0.50 

100 

10 ^ 

43.9 

0.68 

3.11 

0.54 

0.166 

0.70 

0.25 

100 

103.5 

79.2 

0.65 

9.26 

2.21 

0.244 

0.51 

0.15 

100 

10 ^ 

114.7 

0.57 

30.00 

5.30 

0.177 

0.32 

0.17 

10 

104.5 

136.2 

0.50 

98.25 

9.25 

0.092 

0.27 

0.12 

4 


MS30P_SP 


102.5 

21.2 

0.62 

1.23 

0.09 

0.090 

0.79 

0.67 

100 

10 ^ 

43.9 

0.74 

2.92 

0.71 

0.228 

0.65 

0.43 

100 

103.5 

79.2 

0.72 

9.11 

2.32 

0.257 

0.51 

0.32 

100 

10 ^ 

114.7 

0.66 

26.40 

5.70 

0.218 

0.40 

0.27 

10 

104.5 

136.2 

0.55 

96.50 

14.00 

0.145 

0.27 

0.21 

4 


MS3UQ_SP 


102.5 

21.2 

0.58 

1.34 

0.12 

0.097 

0.86 

0.67 

100 

10 ^ 

43.9 

0.73 

3.32 

0.62 

0.177 

0.66 

0.56 

100 

103.5 

79.2 

0.69 

10.28 

2.23 

0.209 

0.59 

0.32 

100 

10 "^ 

114.7 

0.63 

33.60 

6.30 

0.188 

0.45 

0.29 

10 

104.5 

136.2 

0.54 

103.00 

10.50 

0.104 

0.36 

0.26 

4 


MS3S 


102.5 

21.2 

0.44 

1.17 

0.00 

0.000 

0.84 

0.00 

100 

10 ^ 

43.9 

0.55 

2.88 

0.17 

0.043 

0.56 

0.13 

100 

103.5 

79.2 

0.59 

10.40 

1.52 

0.145 

0.22 

0.04 

100 

10 "^ 

114.7 

0.52 

34.30 

3.30 

0.096 

0.06 

0.00 

10 

104.5 

136.2 

0.48 

114.50 

6.25 

0.054 

0.02 

0.04 

4 


MS80P 


102.5 

21.2 

0.93 

1.10 

0.00 

0.000 

0.97 

0.00 

100 

10 ^ 

43.9 

1.02 

2.45 

0.12 

0.036 

0.92 

0.26 

100 

103.5 

79.2 

1.05 

8.81 

0.49 

0.053 

0.54 

0.20 

100 

10 "^ 

114.7 

1.02 

27.90 

1.10 

0.038 

0.39 

0.00 

10 

104.5 

136.2 

1.00 

86.50 

2.25 

0.026 

0.35 

0.00 

4 

105a 

145.3 

0.90 

207.25 

17.50 

0.084 

0.62 

0.05 

4 


MSlOP 


102.5 

21.2 

0.60 

1.16 

0.28 

0.248 

0.91 

0.64 

100 

10^ 

43.9 

0.63 

3.09 

1.48 

0.415 

0.72 

0.35 

100 

103.5 

79.2 

0.49 

9.10 

4.55 

0.516 

0.47 

0.15 

100 

10"^ 

114.7 

0.40 

30.90 

14.30 

0.468 

0.22 

0.05 

10 

104.5 

136.2 

0.33 

106.00 

33.75 

0.320 

0.12 

0.04 

4 

Note. — Columns 1 and 2 show the initial stellar mass 

of the model clusters and the initial maximum stellar mass. 


respectively. Columns 3 and 4 present the averaged half-mass radius,(rh), and the average number of O systems present 
in the calculations, at 3 Myr, respectively. Columns 5 and 6 are the average number of ejected 0-star systems, 

(com and the average 0-star system ejection fraction, {^°’^/ej,o)- The average binary fraction of 0-star systems 

remaining in the cluster, {/b,o)’ average binary fraction of ejected 0-star systems, (/b,ej,o}’ ^ listed 

in columns 7 and 8, respectively. Note that the fraction of binaries among all 0-star systems is smaller than the initial 
value of /b,o = 1 because systems are disrupted in these models. The number of realisations, N’run, for each cluster 
mass is listed in the last column. Only clusters with Med ^ 10^'^ Mq are listed as less massive clusters do not have O 
stars in our models. 

^ The 10^ Mq cluster model is adopted from Banerjee et al. (2012). The initial conditions of the model from Banerjee 
et al. (2012) are slightly different from ours. Their 0-star binaries are mostly close binaries (see Figure 1) and this results 
in higher ejection fractions and less dynamical disruptions of 0-star binaries (higher binary fraction of 0-star systems at 
3 Myr) compared to our models. 



